*** Set working directory ***
cd "C:\..."

*********************** Load Data **********************************************
clear
graph drop _all
cap drop all

import excel "Data.xlsx", sheet("data") firstrow

************************** Setup ***********************************************
* Choose impulse response horizon
local hmax = 36

* Time variable
gen monthly_date = mofd(date)
format monthly_date %tm
drop date
rename monthly_date date
sort date
tsset date
egen time = group(date)
drop if date == .

* Standardize MP Shock series - unit standard deviation
egen mps_std = sd(chg_sveny02)
gen mps = chg_sveny02/mps_std

sort date
tsset date

* Normalize GSCPI for 1998-2019 sample
rename gscpi gscpi_unscl
egen gscpi_mean = mean(gscpi_unscl)
gen gscpi_dm = gscpi_unscl - gscpi_mean
egen gscpi_std = sd(gscpi_dm)
gen gscpi = gscpi_dm/gscpi_std

* GSCPI mp interaction terms
gen gscpi_mps = L.gscpi*mps
gen gscpi_mps1 = L.gscpi*L.mps
gen gscpi_mps2 = L.gscpi*L2.mps

* Put variables in log form
gen cpi_inf = ln(cpi)
gen rgdp_growth = ln(ip)
rename unrt unrt_chg
rename retail retail2
gen retail = ln(retail2)

* Standardize Disaster Shock for 1998-2019 sample
rename costdis costdis_unscl
egen costdis_mean = mean(costdis_unscl)
gen costdis_dm = costdis_unscl - costdis_mean
egen costdis_std = sd(costdis_dm)
gen costdis = costdis_dm/costdis_std

sort date
tsset date

* Interact RHS variables
gen ip_mps = rgdp_growth*mps
gen cpi_mps = cpi_inf*mps
gen unrt_mps = unrt_chg*mps
gen retail_mps = retail*mps
gen ebp_mps = ebp*mps
gen costdis_mps = costdis*mps

* LHS variables
forvalues h = 0/`hmax' {
	gen rgdp_growth_`h' = f`h'.rgdp_growth
	gen cpi_inf_`h' = f`h'.cpi_inf
	gen unrt_chg_`h' = f`h'.unrt_chg
	gen retail_`h' = f`h'.retail
}

************************** LP regressions **************************************
* IP
eststo clear
cap drop b1 u1a d1a u1 d1 Years1 Zero1
gen Years1 = _n-1 if _n<=`hmax'
gen Zero1 =  0    if _n<=`hmax'
gen b1=0
gen u1a=0
gen d1a=0
gen u1=0
gen d1=0

forv h = 0/`hmax' {
	* levels
	newey rgdp_growth_`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(1/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg l(0/2).retail l(0/2).ebp l(0/2).costdis l(1/2).costdis_mps, lag(`h')
replace b1 = _b[gscpi_mps]                    if _n == `h'+1
replace u1a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d1a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u1 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d1 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


* CPI 
cap drop b3 u3a d3a u3 d3 Years2 Zero2
gen Years2 = _n-1 if _n<=`hmax'
gen Zero2 =  0    if _n<=`hmax'
gen b3=0
gen u3a=0
gen d3a=0
gen u3=0
gen d3=0

forv h = 0/`hmax' {
	* levels
	newey cpi_inf_`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(1/2).cpi_inf l(0/2).unrt_chg  l(0/2).retail l(0/2).ebp l(0/2).costdis l(1/2).costdis_mps, lag(`h')
replace b3 = _b[gscpi_mps]                    if _n == `h'+1
replace u3a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d3a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u3 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d3 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


* Unemployment 
cap drop b5 u5a d5a u5 d5 Years3 Zero3
gen Years3 = _n-1 if _n<=`hmax'
gen Zero3 =  0    if _n<=`hmax'
gen b5=0
gen u5a=0
gen d5a=0
gen u5=0
gen d5=0

forv h = 0/`hmax' {
	* levels
	newey unrt_chg_`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(1/2).unrt_chg  l(0/2).retail l(0/2).ebp l(0/2).costdis l(1/2).costdis_mps, lag(`h')
replace b5 = _b[gscpi_mps]                    if _n == `h'+1
replace u5a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d5a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u5 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d5 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


* Retail Sales
cap drop b7 u7a d7a u7 d7 Years4 Zero4
gen Years4 = _n-1 if _n<=`hmax'
gen Zero4 =  0    if _n<=`hmax'
gen b7=0
gen u7a=0
gen d7a=0
gen u7=0
gen d7=0

forv h = 0/`hmax' {
	* levels
	newey retail_`h' l(0/2).mps l.gscpi gscpi_mps gscpi_mps1 gscpi_mps2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg  l(1/2).retail l(0/2).ebp l(0/2).costdis l(1/2).costdis_mps, lag(`h')
replace b7 = _b[gscpi_mps]                    if _n == `h'+1
replace u7a = _b[gscpi_mps] + 1.645* _se[gscpi_mps]  if _n == `h'+1
replace d7a = _b[gscpi_mps] - 1.645* _se[gscpi_mps]  if _n == `h'+1
replace u7 = _b[gscpi_mps] + 1* _se[gscpi_mps]  if _n == `h'+1
replace d7 = _b[gscpi_mps] - 1* _se[gscpi_mps]  if _n == `h'+1

eststo
}


************************* Plot IRF's *******************************************
* IP
twoway ///
(rarea u1a d1a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u1 d1  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b1 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("IP", size(small)) ///

gr rename fig_ip_int, replace


* CPI 
twoway ///
(rarea u3a d3a  Years2,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u3 d3  Years2,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b3 Years2, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero2 Years2, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("CPI", size(small)) ///

gr rename fig_cpi_int, replace


* Unemployment 
twoway ///
(rarea u5a d5a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u5 d5  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b5 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("Unemployment", size(small)) ///

gr rename fig_ur_int, replace


* Retail Sales
twoway ///
(rarea u7a d7a  Years4,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u7 d7  Years4,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b7 Years4, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero4 Years4, lcolor(black)), legend(off) ///
title("Average Supply Chain Conditions", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36, nogrid) ylab(, glpattern(solid)) ///
graphregion(color(white)) plotregion(color(white)) ///
title("Retail Sales", size(small)) ///

gr rename fig_sales_int, replace


* Combine all 4 IRFs
gr combine fig_ip_int fig_cpi_int fig_ur_int fig_sales_int, ///
rows(2) cols(2) graphregion(color(white)) plotregion(color(white)) ///

********************************************************************************
